Sensitivity of rocky planet structures to the equation of state 
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Structures were calculated for Mercury, Venus, Earth, the Moon, and Mars, using a core-mantle 
model and adjusting the core radius to reproduce the observed mass and diameter of each body. 
Structures were calculated using Fe and basalt equations of state of different degrees of sophis- 
tication for the core and mantle. The choice of equation of state had a significant effect on the 
inferred structure. For each structure, the moment of inertia ratio was calculated and compared 
with observed values. Linear Griineisen equations of state fitted to limited portions of shock data 
reproduced the observed moments of inertia significantly better than did more detailed equations 
of state incorporating phase transitions, presumably reflecting the actual compositions of the bod- 
ies. The linear Griineisen equations of state and corresponding structures seem however to be a 
reasonable starting point for comparative simulations of large-scale astrophysical impacts. 



I. INTRODUCTION 

The pressure-temperature-compression equation of 
state (EOS) of condensed matter is vital in understand- 
ing planetary structures, and the response of astrophys- 
ical bodies to impacts. Even for Earth, it is necessary 
to infer the internal structure of the planet from limited 
data. The situation is much more difficult for other bod- 
ies, where seismic data is at best extremely limited and 
more usually absent. 

For planetary structures, we are usually interested in 
pressures from zero to a few hundred gigapascals. In 
this range, EOS were originally deduced empirically from 
shock wave measurements Shock experiments typ- 
ically explore states along the principal Hugoniot of a 
material, which may be significantly hotter than the cor- 
responding compression or pressure occurring in a self- 
gravitating body, and corrections are made to estimate 
the pressure away from the Hugoniot [l]. More recently, 
mechanical presses such as diamond anvil cells have pro- 
vided a source for isothermal compression in this regime, 
although the pressure calibration is made with respect 
to shock measurements on reference materials 0. Sam- 
ples may also be heated within presses, allowing a wide 
range of states to be explored. In parallel with the evo- 
lution of experimental methods, theoretical techniques 
have been developed to predict the EOS from first prin- 
ciples, typically with an a priori accuracy of a few percent 
in mass density, and better if adjusted to reproduce the 
STP state [3[. Theoretical EOS are particularly useful as 
thermodynamically complete states can be calculated for 
arbitrary loading conditions that may be experimentally 
difficult either to induce or to measure. Theoretical EOS 
can also be constructed fairly readily for different mate- 
rial compositions [H . The prediction of phase diagrams is 
however typically less certain than the EOS for a single, 
known phase. 

The strength, failure, and plastic flow behaviors of ma- 
terials are also important in planetary physics. The ef- 
fects of strength must be taken into account when de- 
ducing the EOS from dynamic loading measurements, as 



these are at high strain rates and non-hydrostatic condi- 
tions in contrast to the quasistatic conditions prevailing 
within self-gravitating bodies. Conversely, these proper- 
ties are important for understanding dynamical processes 
including the formation of the body, the effect of impacts, 
cooling, volcanism, mantle convection, and plate tecton- 
ics. These constitutive properties are very challenging 
to predict theoretically, and models rely heavily on ex- 
perimental measurements. Recent studies have demon- 
strated the importance of strength in the heating induced 
by shocks and other dynamic loading Strength and 
flow stress can also change greatly with pressure and tem- 
perature. Recent advances in the measurement and mod- 
eling of plastic flow under extreme conditions has impli- 
cations for astrophysical impacts, such as the interpreta- 
tion of previous simulation studies using simpler mate- 
rial models. An example is the widely-quoted study that 
a liquid interior allow impact energy to be transported 
more efficiently to the antipode of an impact than would 
solid components Q. This conclusion relies strongly on 
the assumptions used for material strength. If the flow 
stress increases sufficiently with compression, solid com- 
ponents can transport compression waves more efficiently 
than can liquids. 

We are constructing theoretical EOS for planetary ma- 
terials of different composition, for use in simulations of 
astrophysical impacts, focusing on systematic differences 
in the shape and location of impact-produced structures 
between the rocky bodies. As a baseline for trial simula- 
tions, we report here the construction of structure mod- 
els for the rocky planets using existing EOS of differ- 
ent sophistication. There is significant uncertainty over 
the composition and temperature profile of all planetary 
bodies; here we use EOS for representative substances 
rather than for the current best estimates of composi- 
tion, because a wider variety of EOS of different degrees 
of sophistication exist for our chosen materials. We op- 
timize a single structure parameter - the core radius - 
for each combination of EOS. This approach is different 
from the commonly-used solution of the Emden equation 
0, [H where the compressibility of the planetary material 
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is assumed to vary linearly with pressure, and the param- 
eters in this relation - i.e., the material properties - are 
deduced from the observable properties of the planet. 



II. ISOSTATIC STRUCTURES 

Consider a spherical body comprising compressible 
material, of mass density p(r). Knowing or assuming 
in addition the distribution of composition and tempera- 
ture, p(r) implies a pressure distribution p(r). The con- 
dition for isostatic equilibrium is that the stress induced 
by pressure variations is balanced by the gravitational 
acceleration g(r) and centripetal force: 



gradp(r) = uJ 2 — p(r)g(r) 



(1) 



where uj is the angular speed. 

For Newtonian gravitation, by Gauss' theorem, g(r) 
can be expressed in terms of the mass enclosed within a 
given radius m(r): 
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m(r) can be calculated simply from the distribution of 
mass density, 



m(r) = An r' 2 p(r') dr' 
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(5) 



In the isostatic equilibrium structure, Eq. [5] can be 
integrated from either the center or the surface r = R to 
obtain the complete structure. In practice, the material 
state at the center depends on the material models used, 
so it is more useful to integrate from the surface, where 
p = and p can be worked out simply for the outermost 
material. 

Usually, the total mass M and outer radius R of the 
body are known fairly well from observation of orbits 
and size respectively, and the problem is to establish the 
internal structure. Starting with the assumption of a 
layered structure, i.e. an ordering of composition, and 
the p{p) relation for each layer, any unknown parameter 
in the structure can be fitted to give the desired total 
mass. 

For example, consider a layered structure where the ra- 
dius of one layer is unknown. For a stable configuration, 
inner layers have a higher density. Consider integration 
downward from the surface as above, assuming a trial 
value for the unknown radius. If too large, the total mass 



M will be reached at some r > 0; if too small, m(0) > 0. 
These two conditions can be used for a solution of this 
shooting problem by bracketing followed by bisection. 

Another remotely observable structural parameter is 
the moment of inertia I, which can be inferred from the 
rate of precession of the rotational pole about the normal 
to the orbit [9(. In the absence of seismic data, I is 
a useful constraint on the structure of the body. The 
moment of inertia for a sphere is 



I = — J p(r)r 4 dr, 



(6) 



and for self-gravitating bodies is conveniently expressed 
as the dimensionless moment of inertia ratio, I /MR 2 . 

Similar approaches assuming the EOS and deducing 
isostatic planetary structures have been reported previ- 
ously. Kerley has calculated structures for Jupiter and 
Saturn, including rotational flattening, using the theory 
of figures and very detailed EOS [IcJ. Seager et al have 
predicted mass-radius relationships for solid exoplanets, 
using EOS which are less suited to subsequent impact 
studies and a different solution approach involvin g in - 
tegration outward from the center of the planet [111 ]. 
Grasset et al have similarly predicted mass-radius rela- 
tionships for rocky and icy exoplanets, again using EOS 
which are less suited to impact studies |12j |. 



III. COMPRESSIBILITY CURVES FOR 
DIFFERENT EQUATIONS OF STATE 

Two types of EOS were considered: linear Griineisen 
EOS using a representative part of the principal shock 
Hugoniot as reference, and more elaborate EOS construc- 
tions incorporating phase transitions and reproducing de- 
tailed features of the shock Hugoniot. The shock Hugo- 
niot data used to calibrate the Griineisen EOS are reliable 
in the sense of being absolute, mechanical measurements, 
but require the additional assumption of a thermal cor- 
rection model to calculate states off the principal Hugo- 
niot, as are required here. 

There is considerable uncertainty in the composition 
of planetary interiors, and detailed EOS are not avail- 
able for all compositions of interest. In order to calculate 
structures for all the rocky planets without introducing 
additional degrees of freedom in the individual compo- 
sitions, the compositions were taken to be pure Fe in 
the core, and basalt in the mantle. At zero pressure, 
basalt typically comprises ~70% plagioclase (often an- 
desite) and ~30% olivine, and is thus reasonably repre- 
sentative of the composition thought to occur through 
the mantle. Further subdivision of layers was not consid- 
ered; the structure chosen allowed the core radius to be 
used as a parameter to ensure the correct total mass for 
each body. Similar structure calculations using detailed 
EOS for other compositions will be reported separately. 

For both Fe and basalt, the shock speed-particle speed 
relation does not follow a straight line over the full pres- 



3 



TABLE I: Linear Griineisen equation of state parameters. 





Fe 


basalt 


Po (g/cm iS ) 


7.86 


2.67 


c (km/s) 


3.635 


1.45 


81 


1.802 


1.97 


r 


2.604 


1.5 



sure range of interest. For basalt, a straight line was cho- 
sen that reproduced the principal Hugoniot at pressures 
of a few tens of gigapascals [l3[ , as a representative range 
for planetary mantles. Similarly, the fit to Fe shock data 
was for pressures ~50-500GPa [13j. The more detailed 
EOS for basalt was a tabulated Griineisen form using 
three piecewise linear fits to the principal Hugoniot [l4l |. 
The detailed Fe EOS was an equilibrium construction 
comprising four solid phases and the fluid, with mixed- 
phase regions [15j . and should be much more accurate 
away from the principal Hugoniot. 

In the absence of accurately-known composition pro- 
files, the use of simple EOS is desirable as it allows plane- 
tary structures to be optimized via the EOS parameters, 
though this optimization was not done here. It is much 
more cumbersome to optimize detailed EOS in this was, 
particularly when presented in tabular form. 

The compression curve used was the principal isen- 
trope. This is a reasonable approach for large, self- 
gravitating objects. One would expect p(p) to be close to 
an isentrope on initial formation, as accretion proceeds 
by multiple small impacts followed by some stratification 
as dense components sink toward the center of the body. 
After formation, the decay of radionuclides (and possibly 
other processes such as tidal deformation and ohmic heat- 
ing) may supply additional heat, but the surface of the 
body cools by radiation. In the adiabatic compression of 
non-porous condensed matter, the thermal contribution 
to pressure is usually much smaller than the repulsion 
between atoms, so the difference between the isentrope 
and the isotherm as extreme limits is small compared 
with the other uncertainties in planetary structure such 
as the composition. 

Each EOS was used in the form p(p, e), and integrated 
numerically starting at STP conditions to obtain p{p) 
along the principal isentrope. The use of the principal 
isentrope for each constituent gives a discontinuity in 
temperature where different layers meet, but the effect on 
pressure is small compared with the uncertainties in com- 
position and temperature. It would be straightforward, 
though less efficient, to calculate the isentrope starting at 
the temperature outside each layer. The isentropes for Fe 
agreed quite closely up to ~120 GPa; at higher pressures, 
the linear Griineisen EOS was significantly stiffer. The 
linear Griineisen EOS for basalt was much softer than 
the tabular EOS at low pressures, and much stiffer at 
high pressures, with a crossover at 100 GPa. (Fig. [T]) 
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FIG. 1: Comparison between isentropes. 



IV. PLANETARY STRUCTURES 

The planetary structure equation, Eq. [5J was inte- 
grated numerically using the 4th order Runge-Kutta al- 
gorithm with a constant step size in radius. The effect 
of rotation was ignored as it depends on latitude; rota- 
tion found to be a negligible effect in a trial equatorial 
integration. The core radius was adjusted to give a den- 
sity profile consistent with the specified total mass of the 
body. The step size was chosen such that core radius was 
converged numerically to 1 km. The resulting moment of 
inertia ratio was converged to 0.1%, which is better than 
measured values. The density near the center of the body 
was very sensitive to small variations in the core radius 
or any other parameters, often giving unphysical varia- 
tions for the innermost couple of tens of kilometers. The 
volume and mass affected were negligible in comparison 
with the overall body. 

Using the linear Griineisen EOS for core and mantle, 
the numerical scheme was very robust, and the core ra- 
dius converged reliably for very wide ranges of the ini- 
tial bracket. Using the tabular EOS, the initial bracket 
had to be chosen much more carefully for a physical core 
radius to be found, because of the increased compress- 
ibility of the core, compressibility of core. Converged 
solutions were however found for all the bodies consid- 
ered and compared with observations of the moment of 
inertia ratio [161 ] and the usually-quoted values of the 
core radius (Table ITT|) . We also summarize pressures cal- 
culated at the key locations of the mantle-core boundary 
and the center (Table |LU| as an indication of the regimes 
desirable to explore experimentally, for each body. 

Given that the EOS used were not for the most likely 
individual core compositions that have been postulated 
(Fe-Ni, or Fe-S for Mars), and the planetary structures 
were simplified, it is interesting to note that the core ra- 
dius inferred for Earth - the only body for which it can be 



TABLE II: Structure parameters for rocky bodies. 



body 


reference 


linear 


tabular 








Griineisen 








y» 


Ij MR 1 


V 


I /MR' 


V 


I/MR Z 




(km) 




(km) 




(km) 




Mercury 


1600 


0.33 


1740 


0.3196 


1802 


0.3140 


Venus 


3000 


0.33 


3014 


0.3348 


3167 


0.3146 


Earth 


3481 


0.33 


3361 


0.3290 


3405 


0.3118 


Moon 


350 


0.393 


620 


0.3647 


752 


0.3572 


Mars 


1700±500 


0.366 


1188 


0.3621 


1706 


0.3316 



TABLE III: Key pressures in rocky bodies. 



body 


linear Griineisen 


tabular 




center 


boundary 


center 


boundary 




(GPa) 


(GPa) 


(GPa) 


(GPa) 


Mercury 


46.5 


9.6 


47.5 


9.6 


Venus 


318 


121 


365 


120 


Earth 


390 


137 


448 


143 


Moon 


8.8 


5.0 


10.0 


4.4 


Mars 


45.2 


27.2 


60.0 


20.7 



estimated using multiple techniques including seismology 
- is too small by only 2% for the tabular EOS and 3% for 
the linear Griineisen. The core radius for the other bod- 
ies did not agree very well with the nominal values in the 
literature, with the exception of Venus, but the difference 
was commensurate with the uncertainty in radius. The 
exceptions were Venus, where the agreement was good at 
0.5-2% (the linear Griineisen EOS being the better), and 
the Moon, where the present calculations were larger by 
a factor ~2. 

For the linear Griineisen EOS, the moment of inertia 
ratios were within 3% except for the Moon, for which it 
was too small by 7%. For the tabular EOS, the ratio 
was systematically smaller: 4-5% except for the Moon 
and Mars (9 and 10% respectively). Too small a ratio 
implies that the body is stiffer than the EOS used in the 
calculation, in the sense that p would be greater for a 
given p. Such an increase in stiffness could reflect a dif- 
ference in composition rather than too low a stiffness for 
the material used. Thus the discrepancy for Mars could 
be attributed to a significant proportion of S, reducing 
the mass density for a given compression. 

The linear Griineisen EOS gives pressures at the center 
of the Earth that are within the range of values in the 
literature (360-390 GPa). 

We next compare in more detail the structures pre- 
dicted for each body by the different EOS. 
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A. Mercury 



For Mercury, the EOS gave very similar pressure dis- 
tributions. The core radius from each EOS was simi- 
lar. The density distribution differed somewhat, partic- 
ularly in the mantle. The pressure to induce the a-e 
phase change occurred just within the core, and showed 
up clearly in the density distribution. The gravitational 
acceleration from each EOS was similar. At the level of 
detail considered here, the EOS were equivalent. (Figs [21 
to SI) 




500 1000 1500 2000 2500 

radius (km) 

FIG. 2: Radial pressure distribution in Mercury calculated 
using different equations of state. 




radius (km) 

FIG. 3: Radial density distribution in Mercury calculated 
using different equations of state. 




radius (km) 



FIG. 4: Radial gravity distribution in Mercury calculated us- 
ing different equations of state. 
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B. Venus 

For Venus, the EOS gave very similar pressure distri- 
butions in the mantle, and deviated monotonically and 
significantly in the core. The core radius from each EOS 
was similar. The density distribution differed somewhat. 
The core pressure for both EOS was well above the a- 
e phase transition. The gravitational acceleration from 
each EOS was similar within the core, but differed some- 
what in the mantle. At the level of detail considered here, 
the EOS were not equivalent but conversely not hugely 
different. (Figs [5] to 0) 
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FIG. 5: Radial pressure distribution in Venus calculated using 
different equations of state. 




FIG. 6: Radial density distribution in Venus calculated using 
different equations of state. 
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C. Earth 

The results for Earth were very similar to those for 
Venus, as one might expect for bodies of similar mass 
and volume. The EOS gave similar pressure distribu- 
tions in the mantle, and deviated monotonically and sig- 
nificantly in the core. The core radius from each EOS 
was very similar. The density distribution differed some- 
what. The core pressure for both EOS was well above 
the a-e phase transition. The gravitational acceleration 
from each EOS differed somewhat over the whole profile 
(except for the ends, which are constrained to be equal), 
although the shapes were very similar. At the level of 
detail considered here, the EOS were not equivalent but 
not hugely different. (Figs 151 to HOI) 
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FIG. 10: Radial gravity distribution in Earth calculated using 
different equations of state. 



FIG. 8: Radial pressure distribution in Earth calculated using 
different equations of state. 



linear Grueneisen 




1000 2000 3000 4000 5000 6000 7000 
radius (km) 



FIG. 9: Radial density distribution in Earth calculated using 
different equations of state. 
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D. Moon 

For the Moon, the EOS gave very similar pressure dis- 
tributions in the mantle, and deviated monotonically and 
significantly in the core. The core radii differed signifi- 
cantly. The density distributions were very similar, par- 
ticularly in the core, except for the difference in core ra- 
dius, differed somewhat, particularly in the mantle. The 
core pressures were all below the a-e phase transition. 
The gravitational acceleration from each EOS was very 
similar in the core, but deviated substantially where the 
core radii differed and in the inner mantle. At the level 
of detail considered here, the EOS were fairly equivalent. 
(Figs ED to EH) 
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FIG. 13: Radial gravity distribution in the Moon calculated 
using different equations of state. 



FIG. 11: Radial pressure distribution in the Moon calculated 
using different equations of state. 
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FIG. 12: Radial density distribution in the Moon calculated 
using different equations of state. 
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E. Mars 

For Mars, the EOS gave similar pressure distributions 
in the mantle, and deviated greatly in the core. The core 
radii differed significantly. The density distributions were 
significantly different, but mostly because of the differ- 
ence in core radius. The core pressures were all above the 
a-e phase transition. The gravitational acceleration from 
each EOS was very similar in the core, but deviated sub- 
stantially where the core radii differed and in the man- 
tle. For simulation purposes, the EOS not equivalent. 
(Figs M to M) 
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FIG. 16: Radial gravity distribution in Mars calculated using 
different equations of state. 



FIG. 14: Radial pressure distribution in Mars calculated using 
different equations of state. 




FIG. 15: Radial density distribution in Mars calculated using 
different equations of state. 
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V. CONCLUSIONS 

The isostatic equations were used to calculate density 
distributions for Mercury, Venus, Earth, the Moon, and 
Mars, assuming different EOS for the structural compo- 
nents, and optimizing the core radius to constrain the 
overall mass to be correct. Two pairs of EOS were con- 
sidered for Fe and basalt, either linear Griineisen fits to 
subsets of published shock data, or more detailed EOS 
constructions including phase changes and nonlinearities 
in the shock data. The core radius of Earth, which is 
by far the best known, was reproduced reasonably well 
using both types of EOS. The other radii, and structures 
for all the bodies, varied significantly between the EOS, 
though to a varying degree. The structures were used to 
calculate the moment of inertia ratio for each body, as 
an independent test of the accuracy of each structure. 

Better overall agreement with the observed proper- 



ties of the rocky bodies considered was obtained with 
the linear Griineisen EOS than the more carefully con- 
structed tabular EOS. This does not imply that the lin- 
ear Griineisen EOS are more accurate models of Fe and 
basalt, but that the properties of the bodies are repre- 
sented more closely by the simpler EOS, presumably be- 
cause the compositions deviate significantly from Fe and 
basalt in the core and mantle respectively. The simpler 
EOS do not reproduce the moment of inertia ratios per- 
fectly, and deviate particularly for the Moon and Mars. 

The linear Griineisen EOS seem to be a reasonable 
starting point for simulations of large-scale impacts on 
these bodies, more so than the tabular EOS. These re- 
sults also suggest that more detailed EOS for the pos- 
tulated compositions of the bodies can be evaluated to 
some degree by their improvement over Fe and basalt in 
calculating the moment of inertia ratios. 
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